% Create Fresh Environment
close all
clear all
clc

L = 201;                % Length of polynomial
t=linspace(-5,5,L);
c = hann(L).*sinc(t)';

% Roots of Unity (how fine to scan unit circle)
mm = 250;                   % nodes in lower half plane
if mm < L/2+1, error('require mm >= N/2 +1')
else M = 2*(mm-1);
end

% Ring precision
rp = 0.03;
rc = 1; %cuurent middle ring
lf_rts=[];
exps = (0:L-1)';            % exponential factors
rr = [rc+rp rc rc-rp];
i = 1;
theta = 0;
indx = 0;
lx = 0;
while(rc > rp)
    
    
    fftA = fft(c.*(rr(1).^exps),M); % evaluate polynomial on the rings
    fftB = fft(c.*(rr(2).^exps),M);
    fftC = fft(c.*(rr(3).^exps),M);

    crA = abs([fftA(end) ; fftA(1:mm+1)]); % append end terms and take abs val
    crB = abs([fftB(end) ; fftB(1:mm+1)]);
    crC = abs([fftC(end) ; fftC(1:mm+1)]);


    k = 2:mm+1;                 % indices of middle-disc without end values
    Bk = crB(k);                % value function on the middle-ring

    % locate the local minima
    locs = (Bk < crA(k-1)) & (Bk < crA(k)) & (Bk < crA(k+1)) ...
                & (Bk < crC(k-1)) & (Bk < crC(k)) & (Bk < crC(k+1)) ...
                    & (Bk < crB(k-1)) & (Bk < crB(k+1));

    indx = find(locs);              % indices of local minima
    lx = length(indx);              % number of minima found

    
    if (lx)
        theta = -2*pi/M*(indx-1);   % convert indices into complex numbers
        lf_rts = [lf_rts; exp(1i*theta)];
    end
    rc = rc - rp;
    rr = [rc+rp rc rc-rp];

end
actual_rts = roots(flipud(c)); % compute roots with companion matrix
plot(lf_rts,'*'), hold on % plot LF approximations,
plot(conj(lf_rts),'*') % companion matrix roots
plot(actual_rts,'or') % and the three concentric discs
uc = exp(1i*linspace(0,2*pi,200)');
plot(uc*rr,'k'), hold off
axis equal, axis off